An External Circular Crack in an Infinite Solid under Axisymmetric Heat Flux Loading in the Framework of Fractional Thermoelasticity

In a real solid there are different types of defects. During sudden cooling, near cracks, there can appear high thermal stresses. In this paper, the time-fractional heat conduction equation is studied in an infinite space with an external circular crack with the interior radius R in the case of axial symmetry. The surfaces of a crack are exposed to the constant heat flux loading in a circular ring R<r<ρ. The stress intensity factor is calculated as a function of the order of time-derivative, time, and the size of a circular ring and is presented graphically.

The classical thermoelasticity is based on the conventional Fourier law for the heat flux vector and the standard parabolic heat conduction equation for temperature. Many theoretical and experimental studies on heat conduction show that in solids with a complex internal structure, in particular with different types of defects, the Fourier law and the parabolic heat conduction equation should be extended to more general relationships. Different generalizations of the Fourier law and heat conduction equation have been studied intensively in the literature (see [24][25][26][27][28][29][30][31] and references therein). The constitutive equations with memory have considerable promise in this area. It should be emphasized that the concept of memory has found wide use in physics, mechanics, economics, and other fields [32][33][34][35][36][37][38][39][40][41][42]. The generalized Fourier law for the heat flux vector q involving memory effects can be written as and leads to the corresponding heat conduction equation with memory, where k and a can be treated as the generalized thermal conductivity and generalized thermal diffusivity coefficient, respectively. The suitable selection of the memory kernel K(t − τ) allows us to obtain different generalized theories of heat conduction and associated generalized thermoelasticity theories. In the case of "full memory" [26,43] , there is no fading of memory, the kernel K(t − τ) is constant, and Equation (2) turns into the wave equation which leads to thermoelasticity without energy dissipation of Green and Naghdi [43]. The exponential memory kernel describes "short-tail memory" and results in the telegraph equation for temperature [24] and generalized thermoelasticity of Lord and Shulman [44] and Green and Lindsay [45].
The theory of integrals and derivatives of fractional order [46][47][48] makes wide use in various areas of science [35,36,[39][40][41][42][49][50][51][52][53][54]. The generalization of the Fourier law with the power "long-tail memory" kernel K(t − τ) [39,[55][56][57] can be expressed in terms of derivatives and integrals of the fractional order q(x, t) = −kI α−1 grad T(x, t), and gives rise to the time-fractional heat conduction equation where [46][47][48] is the Riemann-Liouville fractional integral, denotes the Riemann-Liouville fractional derivative, is the Caputo fractional derivative. Here, Γ(α) is the gamma function. It should be emphasized that in the constitutive equations for the heat flux (3) and (4), the generalized thermal conductivity k has the physical dimension whereas the generalized thermal diffusivity a in the time-fractional heat conduction Equation (7) has the physical dimension This will be important in the following when introducing dimensionless quantities (see Equation (35)).
Theory of thermoelasticity associated with the fractional heat conduction Equation (7) was proposed in [55] (see also the review articles [58,59], and the book [39] which sums up investigations in the field of fractional thermoelasticity). The "middle-tail memory" kernel K(t − τ) in the constitutive Equation (1) is expressed in terms of the Mittag-Leffler functions being the generalization of the exponential function. In this case, we arrive at the time-fractional telegraph equations and the associated theories of thermal stresses [60,61].
Cracks in the framework of generalized theories of thermal stresses were studied in [62][63][64][65][66][67]. In the framework of fractional thermoelasticity, line cracks in a plane [68][69][70] and a penny-shaped crack [71] were considered. In the present paper, we expand the previous studies [68][69][70][71] on the case of an external circular crack with the interior radius R in an infinite solid under the prescribed heat flux at its surfaces. The surfaces of a crack are exposed to the constant heat flux loading in a circular ring R < r < ρ. The stress intensity factor is calculated as a function of the order of time-derivative, time, and the size of a circular ring and is presented graphically.

Formulation of the Problem
We consider a space weakened by an external circular crack with the interior radius R placed in the plane z = 0. The axisymmetric time-fractional heat conduction Equation (7) in the cylindrical coordinates is studied under zero initial conditions and under the prescribed heat flux at the crack surfaces In the subsequent text, we will consider the constant heat flux q 0 acting in the local domain R < r < ρ at the crack surface.
Owing to the geometrical symmetry of the equations with respect to the plane z = 0, we can simplify the problem. Hence, the time-fractional heat conduction equation is investigated in the upper half-space z > 0 under zero initial conditions and under the boundary conditions The zero conditions at infinity are also imposed: In what follows, the Laplace transform with respect to time t will be marked by an asterisk, the Hankel transform of the order zero with respect to the radial coordinate r will be specified by a hat, and the Fourier cosine transform with respect to the coordinate z will be denoted by a tilde (s, ξ, and η are the Laplace, Hankel, and Fourier transform variables, respectively).
The Laplace transform rules for the fractional integrals and derivatives have the following form [46][47][48]:

The Temperature Field
Assuming zero initial conditions for the heat flux, the integral transforms technique gives Here and in the subsequent text, J n (r) denotes the Bessel function of the first kind of the order n, and the following formula for the Fourier cosine transform of the second derivative has been used.
Taking into account that [47,48] where E α,β (z) is the Mittag-Leffler function in two parameters α and β from Equation (29), we obtain Inverting the Hankel and Fourier transforms allows us to obtain the desired solution of the problem (20)-(25) The Mittag-Leffler function E 1,2 (−x) is expressed as [47,53] and after evaluation of the necessary integrals (see Equations (A1) and (A2) in Appendix A), from Equation (34), we obtain the particular case of the solution for the standard parabolic heat conduction equation (α = 1): Figure 1 shows the dependence of temperature on the radial coordinate r for z = 0; the curves presented in Figure 2 describe the dependence of temperature on the spatial coordinate z. In numerical calculations, the following nondimensional quantities have been used (see also Equations (11) and (12)). In the case of the standard heat conduction equation (α = 1), the dimensionless timet reduces to the Fourier number Fo. To simplify numerical calculations, it is convenient to pass to the polar coordinate system in the (ξ,η)-plane:ξ = ζ cos θ,η = ζ sin θ.

The Stress Intensity Factor
The theory of thermoelasticity associated with the fractional heat conduction Equation (7) was proposed in [55]. The system of basic equations of this theory consists of the equilibrium equation in terms of displacements the stress-strain-temperature constitutive equation the geometrical relations and the time-fractional heat conduction equation Here, σ is the stress tensor, e denotes the linear strain tensor, u is the displacement vector, ∇ is the gradient operator, T is the temperature, λ and µ are Lamé constants, K T = λ + 2µ/3 is the bulk modulus, the thermal coefficient of volumetric expansion is designated by β T , I stands for the unit tensor.
As the problem is symmetric with respect to the plane z = 0, we consider the upper half-space z ≥ 0 under the following boundary condition: by virtue of the fact that the surfaces of the external crack are free of mechanical loading, whereas the geometrical symmetry requirements demand that The influence of the temperature field on the stress field can be represented by the displacement potential Φ which is introduced similarly to the classical thermoelasticity [72,73]: with Here, ν denotes the Poisson ratio.
In the axisymmetric case in cylindrical coordinates, the displacement potential gives the components of the displacement vector and the components of the stress tensor where Assuming that ∂Φ(r,z,t) ∂z z=0 = 0, from Equations (33) and (57), we obtain the expression for the displacement potential in the transform domain and after inverting the Hankel and Fourier transforms, we obtain In what follows, we restrict ourselves to calculation of the stress intensity factor which is the most important characteristic of brittle fracture. For the axisymmetric external crack problem in cylindrical coordinates, the stress intensity factor K I describes the singularity of the total stress component σ zz at the crack tip [19] and in terms of the stress component, σ zz is expressed as where From Equations (55) and (59), it is inferred that (63) and Taking into account the integral (A3) from the Appendix A, we arrive at For α = 1, using integral (A4) from the Appendix A, we obtain the particular case of Equation (65) corresponding to the classical thermoelasticity: After passing to the polar coordinate system in the (ξ,η)-plane, Equation (65) is rewritten in nondimensional form as The nondimensional stress factor is presented in Figure 3 as a function of nondimensional timet, and in Figure 4 as a function of the nondimensional parameterρ.
To calculate the total stress field satisfying the boundary conditions (44)- (47), the biharmonic Love function can be used. The interested reader is referred to [73] for general equations and to the paper [71], where this approach was used for a penny-shaped crack.

Concluding Remarks
We have solved the time-fractional heat conduction equation in an infinite solid with an external circular crack with the interior radius R under constant axisymmetric heat flux acting in a circular domain R < r < ρ. The heat flux vector is expressed in terms of the Riemann-Liouville fractional integrals and derivatives of temperature gradient which results in the time-fractional heat conduction equation with the Caputo fractional derivative.
It is worth noting that equations containing the Caputo derivative can be recast to the Riemann-Liouville version (and vice versa) according to the following equation [48,74]: The temperature field and the stress intensity factor are given as integrals with integrands being the Mittag-Leffler functions in two parameters. To evaluate the Mittag-Leffler function E α,2 , we have used the algorithm in [74]. The Mittag-Leffler function E 1/2,2 (−x) is evaluated according to the relation [53] The solution of the classical thermoelasticity problem for the external circular crack is obtained as a particular case when the order of time-derivative α = 1.
For the values 1 < α < 2, we have damped oscillations, but for 2 < α < 3 there arise amplified oscillations. For this reason, the time-fractional heat conduction Equation (7) is considered for the order of fractional derivatives 0 < α ≤ 2.
The time-fractional heat conduction equation for 1 ≤ α ≤ 2 interpolates between the diffusion equation (α = 1) and the wave equation (α = 2) that behave quite differently with respect to their response to a disturbance. The standard heat conduction equation describes a process where a disturbance spreads infinitely fast, whereas in the case of the wave equation the propagation speed of the disturbance is constant. In Figure 1, the bullet points on the r-axis indicate the wave fronts of the solution to the wave equation (α = 2) at r = R − ct and r = ρ + ct, where the coefficient a can be interpreted as a = c 2 with c being the velocity of wave propagation. In terms of dimensionless quantities, the wave fronts arise atr = 1 −t andr =ρ +t.
It is seen from Figure 3 that the stress intensity factor depends significantly on the order of fractional derivatives which corresponds to the power "long-tail memory" kernel K(t − τ). For large values of time, when 0 < α < 1 the stress intensity factor is greater than that for the standard heat conduction (α = 1), whereas for 1 < α ≤ 2 the stress intensity factor is smaller than that for α = 1. The slow heat diffusion (0 < α < 1) increases the stress intensity factor, the fast heat diffusion (1 < α ≤ 2) decreases it. The influence of other memory types, in particular "middle-tail memory"resulting in the time-fractional telegraph equation, will be considered in future research.

Conflicts of Interest:
The authors declare no conflict of interest.